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We use concepts from statistical physics to transform the original evolution equation for the sea 
ice thickness distribution g(h ) due to Thorndike et al. [1] into a Fokker-Planck like conservation 
law. The steady solution is g{h) = Af(q)h q e~ h l H , where q and H are expressible in terms of 
moments over the transition probabilities between thickness categories. The solution exhibits the 
functional form used in observational fits and shows that for h 1, g(h) is controlled by both 
thermodynamics and mechanics, whereas for 1 only mechanics controls g(h). Finally, we derive 
the underlying Langevin equation governing the dynamics of the ice thickness h, from which we 
predict the observed g(h). The genericity of our approach provides a framework for studying the 
geophysical scale structure of the ice pack using methods of broad relevance in statistical mechanics. 
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The Earth’s climate system is a complex nonlinear dy¬ 
namical system [2]. Three main research approaches in 
climate science are common: (a) Observation of the past 
and present state of the system and extrapolation to the 
future, (b) Numerical simulations using Global Circula¬ 
tion Models, which treat the system with the determin¬ 
istic approach of weather forecasting by modeling the 
processes on a coarse-grained scale, (c) Constructing a 
low-order description of the system/subsystem of the cli¬ 
mate, in the vein of theoretical physics. There are sub¬ 
stantial cultural and technical differences between these 
approaches. All have value and all have limitations. Ev¬ 
idently, the ready availability of computing power has 
made approach (c) less favorable. Here we show that one 
of the key variables in polar climate, the sea ice thick¬ 
ness distribution g(h ), can be fruitfully examined quan¬ 
titatively with a core set of tools in statistical physics. 

Although there are hi-fidelity measurements of the area 
of the ice cover during the satellite era [3] , the key quan¬ 
tity reflecting the climatological state of the sea ice cover 
is its volume, making g{h) the central state variable of the 
system. The thickness distribution underlies and reflects 
ice melting or freezing due to the thermodynamic forcing 
of the ocean and atmosphere, and mechanical deforma¬ 
tion; rafting, ridging and the formation of open water 
[1], Nonetheless, although the theory of the ice thickness 
distribution has been with us for forty years we still seek 
a basic understanding of its components in order to test 
its predictions [3]. 

The theory of Thorndike et al. [1] is described by 
a continuous deterministic partial differential equation 
that contains the principal physical processes mentioned 
above and is given by 

f = -V.(u 9 )-|-(/ 9 ) + A (1) 

where u is the velocity of the ice pack and / is its 
growth/melt rate. The principal reason for the diffi¬ 


culty in testing the theory arises from the so-called re¬ 
distribution function, ip, intended to capture mechanical 
processes. Although from observational, theoretical and 
numerical perspectives we have gained a quantitative ex¬ 
planation for many aspects of the the redistribution func¬ 
tion [e.g., 4, 5, and refs, therein], a closed mathematical 
analysis of the original theory is still principally limited 
by this term. In what follows, by viewing ip within the 
framework of kinetic theory, we show that the original 
theory can be rewritten as a Fokker-Planck type equa¬ 
tion. In so doing a number of useful advantages arise. 
First, we can determine the steady solution analytically. 
Second, we provide access to the full range of methods 
and approaches of non-equilibrium statistical physics. In 
this note we describe just two; (1) by comparison with 
observations we deduce the transport coefficients in the 
new evolution equation, which allow for its full numer¬ 
ical solution in a geophysically relevant setting and (2) 
we derive the corresponding Langevin equation for the 
evolution of the ice thickness itself. 

A central ansatz in stochastic dynamics is to consider 
the “microscopic noise” underlying a “macroscopic pro¬ 
cess” as decorrelating on a time scale far faster than the 
macroscopic displacement. The classical test is Brownian 
motion, wherein the inertial macroscopic displacement of 
a pollen grain in water evolves slowly relative to the col- 
lisional events driving its motion. This is embodied in 
the fluctuation-dissipation theorem [e.g., 6]. 

Within this framework we recast ip as follows. Each 
of the “microscopic” mechanical processes that influence 
the ice thickness distribution-rafting, ridging and the for¬ 
mation of open water-occur over a time scale that is very 
rapid relative to the geophysical-scale changes of g{h). 
We thus view these processes as the collisions of solvent 
molecules with a Brownian particle; the individual col¬ 
lisions have a probability of displacing the particle, but 
their phase space is so enormous that we do not study 


2 


them individually. In the same vein, we do not study the 
individual floe-floe interactions in the ice pack. Rather, 
we write ip as 

p OO 

ip = / [g(h + h')w(h + h',h') — g(h)w(h,h')\dh'. (2) 
Jo 

Here, the first and second terms represent the processes 
by which (i) ice floes of thickness h + h' become ice floes 
of thickness h and (ii) ice floes of thickness h become 
ice floes of thickness h — h! respectively, with w dh! the 
transition probability per unit time for these events. 

Taylor expanding the right hand side of equation (2) 
and substituting this is into equation ( 1 ), we obtain 

§ = ~ v • (U9) - §h iM + k (ki 9 ) + w {ki9) - (3) 

where 

poo poo -1 

ki= h'w(h, h!)dh! and k 2 = / —h ,2 w(h,h')dh'. 

Jo Jo ^ 

(4) 

Thus, we have transformed the original theory into a 
Fokker-Planck type of evolution equation. Note that 
in the absence of ice motion equation (3) is exactly the 
Fokker-Planck equation; an advection-diffusion equation 
for the probability density [2, 6 ]. Here, the coefficients 
(Eqs. 4) are the first and second order moments of the 
transition probability between ice thickness categories. 

Choosing L as the horizontal scale, H eq as the vertical 
scale and Uq as the velocity scale, we find three time 
scales in the problem; ( 1 ) the thermal diffusion time scale, 
tjj = H^ q /n, with k the thermal diffusivity of ice; (2) the 
time scale associated with the horizontal motion of ice 
floes, t m = L/Uq ; (3) the relaxation time scale of the 
ice floes when they are involved in collisions, t R ~ 1 / 7 , 
where 7 is the collisional strain rate. These time scales 
are such that, t m ~ t R and r = t R /tD ■C 1. Hence, we 
have f 0 = H eq /t Dl k\ = H eq /t Rl and k 2 = H^ q /t R as the 
scales for the remaining terms. Maintaining the prescaled 
notation, the dimensionless equation can be written in 
one spatial dimension as 

§? = ^ ( " 9 ) - J §h {, 9 )+ th (kl 9 )+ h {k29) ■ (5) 

Now we obtain the steady solely /i-dependent solution 
of equation (5) with boundary conditions g(0) = < 7 ( 00 ) = 
0. The growth rate / in the original theory of Thorndike 
et al. [ 1 ] was determined numerically from the climato- 
logically forced Stefan problem for the ice thickness. If 
AT is the temperature difference over a solid layer of 
thickness h, we take a standard analytical solution for 
its diffusive growth into an isothermal liquid [e.g., 7]. 
Here, one balances heat conduction through the layer 
(oc A T/h) against latent heat production at the inter¬ 
face (oc dh/dt = /) giving 
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with S = L/cpAT the Stefan number, in which L is the 
latent heat of fusion and c p is the specific heat at constant 
pressure. Ignoring the advection term this leads to 



where e = t/S <C 1 because S 1 and r <C 1. Be¬ 
cause the small parameter multiplies regular singulari¬ 
ties, which become 0(1) when h = O(e), we keep all 
terms in equation (7), to which we seek the stationary 
solution and rewrite as 



( 8 ) 


where H = k 2 /k\ and q = e/k 2 . The first integral is 


dg_ 

dh 


+ 



(9) 


where B is the integration constant. We solve equa¬ 
tion (9) using an integrating factor e h / H ~i en ( h ) ^ which 
requires B = 0 to satisfy g(0) = g( 00 ) = 0, and find that 

g(h)=Af(q)h q e- h/H . ( 10 ) 


The prefactor is determined by the normalization condi¬ 
tion fj° g{h)dh = 1 , and is N{q) = [F7 1+9 r(l + g)] \ 
with r(a;) the Euler gamma function. Hence N{q) is 
unique and single valued for M.(q) > —1 and M.(H) > 0. 
Note that q and H have an independent interpretation 
within the framework of the theory and are the sole fit¬ 
ting parameters. Clearly, for h -C 1, g{h) is controlled by 
both thermodynamics and mechanics, whereas for h^$> 1 , 
g{h) is controlled solely by mechanical interactions. 

A Fokker-Planck equation describes the evolution of 
the probability density of a random process, but to study 
the random process itself (h in our case) we study the 
Langevin equation corresponding to equation (7), which 
we write as 

^ = (l-ki) + V^m, (ii) 

where ()) — fci) and y/2k 2 are the drift and diffusion 
terms respectively, and £(i) is Gaussian white noise [e.g., 
8 ]. Clearly, our assumption of k 2 being a constant in 
the determination of the solution for g(h ) in equation 
( 10 ), translates into additive noise in the corresponding 
Langevin equation shown in (11). 

We now compare our theory with the thickness distri¬ 
butions obtained during the ICESat mission [9]. Figure 
1 shows fits of our solution ( 10 ) to the distribution func¬ 
tions for the period February-March (F-M) for (a) 2008 
and (b) 2004, which we have chosen to demonstrate both 
the typical fit (2008) and the worst fit (2004) of our so¬ 
lution to the observations [c.f., Fig. 6 of 9]. The key 
reasons for deviations are ( 1 ) the observations span the 






3 



FIG. 1. Comparison of our theory with satellite measure¬ 
ments for February-March (F-M) of (a) 2008 and (b) 2004. 
Circles are the distribution functions from ICESat [9] and 
lines are the fits using equation (10). In (a), q = 1.849 and 
H = 0.788 m. and in I'M. a = 1.848 and ff = 0.910 m. 




FIG. 2. Solution to the Langevin equation (11) for F-M 
2008. Here, e = 0.046, ki = 0.048 and k 2 = 0.025. The 
ensemble size used is N en = 10 5 and the time step for the 
Euler-Maruyama scheme [10] is At = 10~ 5 . The total non- 
dimensional integration time is T = 140. Figures (a) and 
(b) show the thickness distribution and four realizations from 
the ensemble respectively. In (a), the circles represent g(h) 
from the Langevin equation and the solid curve is the fit using 
equation (10), which gives q = 1.838 and H = 0.777 m. 

ice cover, and yet near landmasses, and depending on 
wind direction, it is possible that k\ and k 2 will differ 
locally, (2) we neglected the advection term when deter¬ 
mining the solution, (3) the form of / used in equation 
(3) is the solution of the ideal Stefan problem for growth 
only, and (4) we are comparing the steady solution to the 
data. Incorporating these and related issues are key as¬ 
pects of a thorough numerical study of equation 5, which 
are part of a longer treatment. 

Now from the values of q and H we obtain k\ and 
k 2 and use them in equation (11) to evolve h itself. The 
solution to the Langevin equation corresponding to figure 
1(a) is shown in figure 2. Invoking ergodicity, clearly 
equation (11) qualitatively reproduces the observed g(h) 
and thus acts as an ideal and simple model to study the 
thickness distribution. 

We have transformed the original evolution equa¬ 
tion for the sea ice thickness distribution, g(h ), due to 
Thorndike et al. [1], to a Fokker-Planck like equation by 
recasting the redistribution function, ip, using the anal¬ 


ogy with the theory Brownian motion. The idea is that 
the mechanical processes embodied in ip (rafting, ridging 
and the formation of open water) are thought of like the 
collisions of solvent molecules with a Brownian particle- 
the individual events that change the ice thickness occur 
on time and space scales that are short relative to the 
geophysical-scale changes of g(h). Thus, we do not treat 
the individual floe-floe interactions in the ice pack, but 
rather only the moments of the transition probabilities 
for these events. That the integrals describing these mo¬ 
ments rapidly converge to take constant values is borne 
out by comparison with observations; the stationary so¬ 
lution (equation 10 ) of the new evolution equation cap¬ 
tures the basin-scale satellite measurements of the dis¬ 
tribution. Finally, the corresponding Langevin equation 
( 11 ) is evolved with observationally constrained param¬ 
eters to study the evolution of h itself. The associated 
agreement of the g(h ) obtained from this approach with 
the observations is consistent with an ergodic thickness 
field. The simplicity of the approach and its immediate 
connection with the edifice of non-equilibrium statistical 
mechanics make it appealing for a wide range of reasons, 
from context for comparison with observations to simpli¬ 
fication of models. 
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